function h = grid_spy(grid,varargin)
% h = grid_spy(grid,dofs,ddc_grid,ddc_dofs)
%
% Plot grid connectivity information as follows:
%  vertices : o, black
%  sides    : x, green, label aligned with side normal
%  elements : blue, a square marks the first local vertex
%  boundary : boundary markers displayed in magenta

 defined_dofs = 0;
 defined_ddc_g = 0;
 defined_ddc_d = 0;
 if(nargin>1)
   dofs = varargin{1};
   if(~isempty(dofs))
     defined_dofs = 1;
   end
   if(nargin>2)
     ddc_grid = varargin{2};
     if(~isempty(ddc_grid))
       defined_ddc_g = 1;
     end
     if(nargin>3)
       ddc_dofs = varargin{3};
       if(~isempty(ddc_dofs))
         defined_ddc_d = 1;
       end
     end
   end
 end

 h = figure;

 delta = 0.03*(grid.vol/factorial(grid.d))^(1/grid.d);
 done = 0;
 while(not(done))
   disp('Insert a string with the required entity:')
   disp('  B: boundary (red)')
   disp('  V: vertexes (blue)')
   disp('  E: elements (black)')
   disp('  L: sides (magenta)')
   disp('  N: nodes (green)')
   disp('  d: neighboring subdomains - vertexes and sides (red)')
   disp('  i: neighboring indexes - vertexes and sides (gray)')
   disp('  v: global vertexes (blue)')
   disp('  l: global sides (magenta)')
   disp('  e: global elements (black)')
   disp('  n: neighboring nodes (green)')
   disp('  t: node neighbours (green)')
   disp('  g: global nodes (green)')
   disp('  Q: quit')
   disp('')
   disp('Examples:')
   disp('  BE   -> boundary and elements')
   disp('  EL   -> elements and sides')
   disp('  BELN -> boundary, elements, sides e nodes')

   fflush(stdout);
   string = input('selection: ');
   if(not(isempty(find(string=='Q'))))
     return
   end

   % grid
   hold off
   x_ed = zeros(grid.d,2,grid.me.children(2).ns,grid.ne);
   for ie=1:grid.ne
     for ic=1:grid.me.children(2).ns
       for ivl=1:2
         iv = grid.e(ie).iv(grid.me.children(2).s(ivl,ic));
         x_ed(:,ivl,ic,ie) = grid.v(iv).x;
       end
     end
   end
   switch grid.d
    case 1
     error('Not implemented')
    case 2
     x_ed_p = reshape(x_ed(1,:,:,:),2,grid.me.children(2).ns*grid.ne);
     y_ed_p = reshape(x_ed(2,:,:,:),2,grid.me.children(2).ns*grid.ne);
     plot(x_ed_p,y_ed_p,'-y');
     xlabel('x'), ylabel('y')
    case 3
     x_ed_p = reshape(x_ed(1,:,:,:),2,grid.me.children(2).ns*grid.ne);
     y_ed_p = reshape(x_ed(2,:,:,:),2,grid.me.children(2).ns*grid.ne);
     z_ed_p = reshape(x_ed(3,:,:,:),2,grid.me.children(2).ns*grid.ne);
     plot3(x_ed_p,y_ed_p,z_ed_p,'-y');
     xlabel('x'), ylabel('y'), zlabel('z')
   end
   axis equal
   hold on

   v_delta = zeros(grid.nv,1);
   s_delta = zeros(grid.ns,1);
   e_delta = zeros(grid.ne,1);
   if(defined_dofs)
     n_delta = zeros(dofs.ndofs,1);
   end
   for i=1:length(string)
     switch(string(i))

      case('B')
       for ib=1:grid.nb
         is  = ib + grid.ni;
	 ie  = grid.s(is).ie(1);
	 isl = grid.s(is).isl(1);
	 pos = grid.s(is).xb + ...
	         0.2*get_h(grid.s(is))*grid.e(ie).n(:,isl);
         switch grid.d
          case 2
           text(pos(1),pos(2),num2str(-grid.s(is).ie(2)), ...
	        'color','red','FontSize',14);
          case 3
           text(pos(1),pos(2),pos(3),num2str(-grid.s(is).ie(2)), ...
	        'color','red','FontSize',14);
	 end
       end

      case('V')
       for iv=1:grid.nv
         pos = grid.v(iv).x;
         pos(2) = pos(2)+v_delta(iv);
         switch grid.d
          case 2
           text(pos(1),pos(2),num2str(iv), ...
	        'color','blue','FontSize',9);
          case 3
           text(pos(1),pos(2),pos(3),num2str(iv), ...
	        'color','blue','FontSize',9);
	 end
         v_delta(iv)  = v_delta(iv) + delta;
       end

      case('E')
       for ie=1:grid.ne
         pos = grid.e(ie).xb;
         pos(2) = pos(2)+e_delta(ie);
         switch grid.d
          case 2
           text(pos(1),pos(2),num2str(ie), ...
	        'color','black','FontSize',9);
          case 3
           text(pos(1),pos(2),pos(3),num2str(ie), ...
	        'color','black','FontSize',9);
	 end
         e_delta(ie) = e_delta(ie) + delta;
       end

      case('L')
       for is=1:grid.ns
         pos = grid.s(is).xb;
         pos(2) = pos(2)+s_delta(is);
         switch grid.d
          case 2
           text(pos(1),pos(2),num2str(is), ...
	        'color','magenta','FontSize',9);
          case 3
           text(pos(1),pos(2),pos(3),num2str(is), ...
	        'color','magenta','FontSize',9);
	 end
         s_delta(is)  = s_delta(is) + delta;
       end

      case('N')
       if(defined_dofs)
         for in=1:dofs.ndofs
           pos = dofs.x(:,in);
           switch grid.d
            case 2
	     plot(pos(1),pos(2),'o','color','green','MarkerSize',1);
             pos(1) = pos(1)+0.1*delta;
             pos(2) = pos(2)+n_delta(in);
             text(pos(1),pos(2),num2str(in), ...
                  'color','green','FontSize',9);
            case 3
	     plot3(pos(1),pos(2),pos(3),'o','color','green','MarkerSize',1);
             pos(1) = pos(1)+0.1*delta;
             pos(2) = pos(2)+n_delta(in);
             pos(3) = pos(3)+0.1*delta;
             text(pos(1),pos(2),pos(3),num2str(in), ...
                  'color','green','FontSize',9);
           end
           n_delta(in)  = n_delta(in) + delta;
         end
       else
         warning('You need to provide the optional input dofs')
       end

      case('d')
       if(defined_ddc_g)
         for ivb=1:ddc_grid.nnv
	   iv = ddc_grid.nv(ivb).i;
	   id = ddc_grid.nv(ivb).id;
           pos = grid.v(iv).x;
           pos(2) = pos(2)+v_delta(iv);
           switch grid.d
            case 2
             text(pos(1),pos(2),num2str(id), ...
                  'color',[0.6 0.2 0.2],'FontSize',9);
            case 3
           text(pos(1),pos(2),pos(3),num2str(id), ...
                  'color',[0.6 0.2 0.2],'FontSize',9);
           end
           v_delta(iv) = v_delta(iv) + delta;
         end
         for isb=1:ddc_grid.nns
	   is = ddc_grid.ns(1,isb);
	   id = ddc_grid.ns(2,isb);
           pos = grid.s(is).xb;
           pos(2) = pos(2)+s_delta(is);
           switch grid.d
            case 2
             text(pos(1),pos(2),num2str(id), ...
                  'color',[0.6 0.2 0.2],'FontSize',9);
            case 3
           text(pos(1),pos(2),pos(3),num2str(id), ...
                  'color',[0.6 0.2 0.2],'FontSize',9);
           end
           s_delta(is)  = s_delta(is) + delta;
         end
       else
         warning('You need to provide the optional input ddc_grid')
       end

      case('i')
       if(defined_ddc_g)
         for ivb=1:ddc_grid.nnv
	   iv = ddc_grid.nv(ivb).i;
	   in = ddc_grid.nv(ivb).in;
           pos = grid.v(iv).x;
           pos(2) = pos(2)+v_delta(iv);
           switch grid.d
            case 2
             text(pos(1),pos(2),num2str(in), ...
                  'color',[0.4 0.4 0.4],'FontSize',9);
            case 3
           text(pos(1),pos(2),pos(3),num2str(in), ...
                  'color',[0.4 0.4 0.4],'FontSize',9);
           end
           v_delta(iv) = v_delta(iv) + delta;
         end
         for isb=1:ddc_grid.nns
	   is = ddc_grid.ns(1,isb);
	   in = ddc_grid.ns(3,isb);
           pos = grid.s(is).xb;
           pos(2) = pos(2)+s_delta(is);
           switch grid.d
            case 2
             text(pos(1),pos(2),num2str(in), ...
                  'color',[0.4 0.4 0.4],'FontSize',9);
            case 3
           text(pos(1),pos(2),pos(3),num2str(in), ...
                  'color',[0.4 0.4 0.4],'FontSize',9);
           end
           s_delta(is)  = s_delta(is) + delta;
         end
       else
         warning('You need to provide the optional input ddc_grid')
       end

      case('v')
       if(defined_ddc_g)
         for iv=1:grid.nv
           pos = grid.v(iv).x;
           pos(2) = pos(2)+v_delta(iv);
	   ig = ddc_grid.gv(1,iv);
           switch grid.d
            case 2
             text(pos(1),pos(2),num2str(ig), ...
                  'color',[0.2 0.2 0.6],'FontSize',9);
            case 3
           text(pos(1),pos(2),pos(3),num2str(ig), ...
                  'color',[0.2 0.2 0.6],'FontSize',9);
           end
           v_delta(iv) = v_delta(iv) + delta;
         end
       else
         warning('You need to provide the optional input ddc_grid')
       end

      case('l')
       if(defined_ddc_g)
         for is=1:grid.ns
           pos = grid.s(is).xb;
           pos(2) = pos(2)+s_delta(is);
	   ig = ddc_grid.gs(1,is);
           switch grid.d
            case 2
             text(pos(1),pos(2),num2str(ig), ...
                  'color',[0.8 0.1 0.6],'FontSize',9);
            case 3
           text(pos(1),pos(2),pos(3),num2str(ig), ...
                  'color',[0.8 0.1 0.6],'FontSize',9);
           end
           s_delta(is) = s_delta(is) + delta;
         end
       else
         warning('You need to provide the optional input ddc_grid')
       end

      case('e')
       if(defined_ddc_g)
         for ie=1:grid.ne
           pos = grid.e(ie).xb;
           pos(2) = pos(2)+e_delta(ie);
	   ig = ddc_grid.ge(ie);
           switch grid.d
            case 2
             text(pos(1),pos(2),num2str(ig), ...
                  'color',[0.2 0.2 0.2],'FontSize',9);
            case 3
           text(pos(1),pos(2),pos(3),num2str(ig), ...
                  'color',[0.2 0.2 0.2],'FontSize',9);
           end
           e_delta(ie) = e_delta(ie) + delta;
         end
       else
         warning('You need to provide the optional input ddc_grid')
       end

      case('n')
       if(defined_ddc_d)
         for ivb=1:ddc_dofs.nndofs
	   iv = ddc_dofs.ndofs(ivb).i;
	   in = ddc_dofs.ndofs(ivb).in;
           pos = dofs.x(:,iv);
           switch grid.d
            case 2
	     plot(pos(1),pos(2),'o', ...
	          'color','green','MarkerSize',1);
             pos(1) = pos(1)+0.1*delta;
             pos(2) = pos(2)+n_delta(iv);
             text(pos(1),pos(2),num2str(in), ...
                  'color',[0.3 0.8 0.2],'FontSize',9);
            case 3
           text(pos(1),pos(2),pos(3),num2str(in), ...
                  'color',[0.3 0.8 0.2],'FontSize',9);
           end
           n_delta(iv) = n_delta(iv) + delta;
         end
       else
         warning('You need to provide the optional input ddc_grid')
       end

      case('t')
       if(defined_ddc_d)
         for ivb=1:ddc_dofs.nndofs
	   iv = ddc_dofs.ndofs(ivb).i;
	   in = ddc_dofs.ndofs(ivb).id;
           pos = dofs.x(:,iv);
           switch grid.d
            case 2
	     plot(pos(1),pos(2),'o', ...
	          'color','green','MarkerSize',1);
             pos(1) = pos(1)+0.1*delta;
             pos(2) = pos(2)+n_delta(iv);
             text(pos(1),pos(2),num2str(in), ...
                  'color',[0.1 0.4 0.1],'FontSize',9);
            case 3
           text(pos(1),pos(2),pos(3),num2str(in), ...
                  'color',[0.1 0.4 0.1],'FontSize',9);
           end
           n_delta(iv) = n_delta(iv) + delta;
         end
       else
         warning('You need to provide the optional input ddc_grid')
       end

      case('g')
       if(defined_ddc_d)
         for in=1:dofs.ndofs
           pos = dofs.x(:,in);
           switch grid.d
            case 2
	     plot(pos(1),pos(2),'o', ...
	          'color','green','MarkerSize',1);
             pos(1) = pos(1)+0.1*delta;
             pos(2) = pos(2)+n_delta(in);
             text(pos(1),pos(2),num2str(ddc_dofs.gdofs(1,in)), ...
                  'color',[0.2 0.7 0.2],'FontSize',9);
            case 3
           text(pos(1),pos(2),pos(3),num2str(ddc_dofs.gdofs(1,in)), ...
                  'color',[0.2 0.7 0.2],'FontSize',9);
           end
           n_delta(in)  = n_delta(in) + delta;
         end
       else
         warning('You need to provide the optional input ddc_dofs')
       end

      otherwise
       warning(['Unknown entity "',string(i),'".'])

     end
   end
 end
return

% my_pdemesh(h,grid,"y");
% hold on
% axis normal
%
% xmin = +inf;
% xmax = -inf;
% for iv=1:grid.nv
%   xmin = min(xmin,min(grid.v(iv).x));
%   xmax = max(xmax,max(grid.v(iv).x));
% end
% dx = (xmax-xmin)/(8*sqrt(grid.ne));
%
% % vertex information
% for iv=1:grid.nv
%  switch(grid.d)
%   case 1
%    hh = plot(grid.v(iv).x(1),0,'ko');
%    set(hh,'markersize',3);
%    text(grid.v(iv).x(1),dx,num2str(iv),'color','k');
%   case 2
%    hh = plot(grid.v(iv).x(1),grid.v(iv).x(2),'ko');
%    set(hh,'markersize',3);
%    text(grid.v(iv).x(1)+dx,grid.v(iv).x(2)+dx,num2str(iv),'color','k');
%   case 3
%    hh = plot3(grid.v(iv).x(1),grid.v(iv).x(2),grid.v(iv).x(3),'ko');
%    set(hh,'markersize',3);
%    text(grid.v(iv).x(1)+dx,grid.v(iv).x(2)+dx,grid.v(iv).x(3)+dx, ...
%      num2str(iv),'color','k');
%  end
% end
%
% % element information
% for ie=1:grid.ne
%  switch(grid.d)
%   case 1
%    text(grid.e(ie).xb(1),dx,num2str(ie),'color','b');
%    x0 = 0.6*grid.e(ie).x0 + 0.4*grid.e(ie).xb;
%    hh = plot(x0(1),0,'bs');
%    set(hh,'markersize',3);
%   case 2
%    text(grid.e(ie).xb(1),grid.e(ie).xb(2),num2str(ie),'color','b');
%    x0 = 0.6*grid.e(ie).x0 + 0.4*grid.e(ie).xb;
%    hh = plot(x0(1),x0(2),'bs');
%    set(hh,'markersize',3);
%   case 3
%    text(grid.e(ie).xb(1),grid.e(ie).xb(2),grid.e(ie).xb(3), ...
%      num2str(ie),'color','b');
%    x0 = 0.6*grid.e(ie).x0 + 0.4*grid.e(ie).xb;
%    hh = plot3(x0(1),x0(2),x0(3),'bs');
%    set(hh,'markersize',3);
%    x0 = 0.6*grid.v(grid.e(ie).iv(2)).x + 0.4*grid.e(ie).xb;
%    hh = plot3(x0(1),x0(2),x0(3),'b*');
%    set(hh,'markersize',3);
%  end
% end
%
% % side information
% for is=1:grid.ns
%  switch(grid.d)
%   case 1
%    hh = plot(grid.s(is).xb(1),0,"gx");
%    set(hh,'markersize',3);
%    hh = text(grid.s(is).xb(1)+0.3*dx*grid.s(is).n(1), ...
%         -dx,num2str(is),'color','g');
%    set(hh,'rotation',180/pi*atan2(0,grid.s(is).n(1)));
%   case 2
%    hh = plot(grid.s(is).xb(1),grid.s(is).xb(2),"gx");
%    set(hh,'markersize',3);
%    hh = text(grid.s(is).xb(1)+0.3*dx*grid.s(is).n(1), ...
%         grid.s(is).xb(2)+0.3*dx*grid.s(is).n(2),num2str(is),'color','g');
%    set(hh,'rotation',180/pi*atan2(grid.s(is).n(2),grid.s(is).n(1)));
%   case 3
%    hh = plot3(grid.s(is).xb(1),grid.s(is).xb(2),grid.s(is).xb(3),"gx");
%    set(hh,'markersize',3);
%    hh = text(grid.s(is).xb(1)+0.3*dx*grid.s(is).n(1), ...
%              grid.s(is).xb(2)+0.3*dx*grid.s(is).n(2), ...
%              grid.s(is).xb(3)+0.3*dx*grid.s(is).n(3), ...
%	      num2str(is),'color','g');
%  end
% end
% % boundary side information
% for is=grid.ni+1:grid.ns
%  switch(grid.d)
%   case 1
%    text(grid.s(is).xb(1)+5*dx*grid.s(is).n(1),0, ...
%	 num2str(-grid.s(is).ie(2)),'color','m','FontSize',15);
%   case 2
%    text(grid.s(is).xb(1)+5*dx*grid.s(is).n(1), ...
%         grid.s(is).xb(2)+5*dx*grid.s(is).n(2), ...
%	 num2str(-grid.s(is).ie(2)),'color','m','FontSize',15);
%   case 3
%    text(grid.s(is).xb(1)+5*dx*grid.s(is).n(1), ...
%         grid.s(is).xb(2)+5*dx*grid.s(is).n(2), ...
%         grid.s(is).xb(3)+5*dx*grid.s(is).n(3), ...
%	 num2str(-grid.s(is).ie(2)),'color','m','FontSize',15);
%  end
% end
%
% axis equal


function h = get_h(s)
 d = s.d;
 if(d==1)
   h = 1;
 else
   h = (factorial(d-1)*sqrt(2^(d-1)/d)*s.a)^(1/(d-1));
 endif
return

